#include "AnalysisPipeline.h"

AnalysisPipeline::AnalysisPipeline(int tSampleNum, float E, float I, int L, int tTaxonomicLevel)
{
    Threshold=0.5;
    SampleNum=tSampleNum; 
    BLASTEvalueThreshold=E;
    BLASTIdentityThreshold=I;
    BLASTMappingLengthThreshold=L;
    ThetaIndexValue=NULL;
    TaxonomicLevel=tTaxonomicLevel;
}

AnalysisPipeline::~AnalysisPipeline()
{
    map<string,double*>::iterator it=Abundance.begin();
    while(it!=Abundance.end())
    {
        delete[] (*it).second;
         it++;
    }
    Abundance.clear();
    it=WilcoxonScore.begin();
    while(it!=WilcoxonScore.end())
    {
        if((*it).second!=NULL)
            delete[] (*it).second;
        it++;
    }
    WilcoxonScore.clear();
    RawMapping.clear();
    TaxonomyID.clear();
    rTaxonomyID.clear();
    QueryID.clear();
    for(int i=0;i<SampleNum;i++)
    {
        delete[] ThetaIndexValue[i];
        delete[] BrayCurtisDissimilarity[i];
        delete[] JaccardIndexValue[i];
        delete[] SensitivityMatrix[i];
    }
    delete[] ThetaIndexValue;
    delete[] BrayCurtisDissimilarity;
    delete[] JaccardIndexValue;
    delete[] SensitivityMatrix;
    delete[] BLASTMapping;
    delete[] DataSize;
    delete[] AmbiSize;
    if(Median!=NULL)
        delete[] Median;
}

int AnalysisPipeline::infer_taxonomy(int OperationFlag, bool MedianNormalizationFlag)
{
    DataSize=new double[SampleNum];
    AmbiSize=new int[SampleNum];
    memset(DataSize,0,sizeof(double)*SampleNum);
    for(int i=0;i<SampleNum;i++)
    {
        AmbiSize[i]=BLASTMapping[i].AmbiCount;
        BLASTMapping[i].Iterator_reset();
        while(!BLASTMapping[i].Iterator_end())
        {
            DataSize[i]+=(*BLASTMapping[i].Iterator).second;
            if(OperationFlag==2)
            {
                map<string,double*>::iterator it = Abundance.find((*BLASTMapping[i].Iterator).first);
                if (it == Abundance.end())
                {
                    Abundance[(*BLASTMapping[i].Iterator).first]=new double[SampleNum];
                    memset(Abundance[(*BLASTMapping[i].Iterator).first],0,sizeof(double)*SampleNum);
                }
                Abundance[(*BLASTMapping[i].Iterator).first][i]+=(*BLASTMapping[i].Iterator).second;
            }
            else
            {
                map<string,double*>::iterator it = Abundance.find(GenBank.get_Taxonomy((*BLASTMapping[i].Iterator).first,TaxonomicLevel));
                if (it == Abundance.end())
                {
                    Abundance[GenBank.get_Taxonomy((*BLASTMapping[i].Iterator).first, TaxonomicLevel)]=new double[SampleNum];
                    memset(Abundance[GenBank.get_Taxonomy((*BLASTMapping[i].Iterator).first, TaxonomicLevel)],0,sizeof(double)*SampleNum);
                }
                Abundance[GenBank.get_Taxonomy((*BLASTMapping[i].Iterator).first, TaxonomicLevel)][i]+=(*BLASTMapping[i].Iterator).second;
            }
            BLASTMapping[i].Iterator++;
        }
    }
    if(MedianNormalizationFlag==false)
    {
        map<string,double*>::iterator it = Abundance.begin();
        while(it!=Abundance.end())
        {   
            for(int i=0;i<SampleNum;i++)
                (*it).second[i]/=DataSize[i];
            it++;
        }
    }
    else
        normalize_median();
    cout<<"Infer Taxonomy Complete"<<endl;
    return 1;
}

void AnalysisPipeline::dist_calc(int SampleIndex)
{
     cout<<"Loading distance matrix\n";
     Dist=new double*[BLASTMapping[SampleIndex].MappingIterator_size()];
     T=new double*[BLASTMapping[SampleIndex].MappingIterator_size()];
     BLASTMapping[SampleIndex].MappingIterator_reset();
     int count=0;
     while(!BLASTMapping[SampleIndex].MappingIterator_end())
     {
         QueryID[(*BLASTMapping[SampleIndex].MappingIterator).first]=count;
         Dist[count]=new double[GenBank.Iterator_size(2)];
         T[count]=new double[GenBank.Iterator_size(2)];
         for(int i=0;i<GenBank.Iterator_size(2);i++)
             Dist[count][i]=0.0;
         memset(T[count],0,sizeof(double)*(GenBank.Iterator_size(2)));
         count++;
         BLASTMapping[SampleIndex].MappingIterator++;
     }
     BLASTMapping[SampleIndex].MappingIterator_reset();
     while(!BLASTMapping[SampleIndex].MappingIterator_end())
     {
         map<string, MappingInfo>::iterator it=(*BLASTMapping[SampleIndex].MappingIterator).second.begin();
         while(it!=(*BLASTMapping[SampleIndex].MappingIterator).second.end())
         {
             Dist[QueryID[(*BLASTMapping[SampleIndex].MappingIterator).first]][TaxonomyID[GenBank.get_Taxonomy((*it).first, TaxonomicLevel)]]=MAX(Dist[QueryID[(*BLASTMapping[SampleIndex].MappingIterator).first]][TaxonomyID[GenBank.get_Taxonomy((*it).first, TaxonomicLevel)]],(double)gsl_ran_gaussian_pdf(((*it).second.Identity/100.0-1.0), Threshold));
             it++;
         }
         BLASTMapping[SampleIndex].MappingIterator++;
     }
     for(int i=0;i<QueryID.size();i++)
     {
         double total=0;
         for(int j=0;j<TaxonomyID.size();j++)
             total+=Dist[i][j];
         for(int j=0;j<TaxonomyID.size();j++)
             Dist[i][j]/=total;
     }
     cout<<"Loading distance matrix done\n";
     return;
}

int AnalysisPipeline::EM(int SampleIndex)
{
    cout<<"Using EM algorithm to infer abundance for sample "<<SampleIndex<<"\n";
    //simple initial
    InitialGuess=new double[GenBank.Iterator_size(2)];
    double* PrevGuess=new double[GenBank.Iterator_size(2)];
    for(int i=0;i<TaxonomyID.size();i++)
        InitialGuess[i]=(double)1.0/TaxonomyID.size();
    double threshold=0;
    int count=0;
    do
    {   //E-step
        threshold=0;
        for(int i=0;i<QueryID.size();i++)
        {
            double sum=0;
            for(int j=0;j<TaxonomyID.size();j++)
            {
                T[i][j]=InitialGuess[j]*Dist[i][j];
                sum+=T[i][j];
            }
            for(int j=0;j<TaxonomyID.size();j++)
                T[i][j]/=sum;
        }
        //M-step
        memcpy(PrevGuess, InitialGuess, sizeof(double)*(TaxonomyID.size()));
        memset(InitialGuess, 0, sizeof(double)*(TaxonomyID.size()));
        for(int i=0;i<TaxonomyID.size();i++)
        {
            for(int j=0;j<QueryID.size();j++)
            {
                InitialGuess[i]+=T[j][i];
            }
            InitialGuess[i]/=(double)QueryID.size();
            threshold=MAX(threshold, (double)abs(InitialGuess[i]-PrevGuess[i])/PrevGuess[i]);
        }
        count++;        
    }while(threshold>0.01&&count<100);
    
    /*
    //final assignment
    memset(InitialGuess, 0, sizeof(double)*(TaxonomyID.size()));
    for(int i=0;i<QueryID.size();i++)
    {
        double MAX=-1;
        int Index=-1;
        for(int j=0;j<TaxonomyID.size();j++)
        {
            if(T[i][j]>MAX)
            {
                MAX=T[i][j];
                Index=j;
            }
        }
        InitialGuess[Index]+=1.0;
    }
    for(int i=0;i<TaxonomyID.size();i++)
        InitialGuess[i]/=(double)QueryID.size();
    */
    
    for(int i=0;i<BLASTMapping[SampleIndex].MappingIterator_size();i++)
    {
        delete[] Dist[i];
        delete[] T[i];
    }
    delete[] Dist;
    delete[] T;
    delete[] PrevGuess;
    QueryID.clear();
    cout<<"EM complete\n";
    return 1;
}

int AnalysisPipeline::infer_taxonomy_EM()
{
     int count=0;
     GenBank.Iterator_reset(2);
     while(!GenBank.Iterator_end(2))
     {
         TaxonomyID[(*GenBank.GBIterator).first]=count;
         rTaxonomyID.push_back((*GenBank.GBIterator).first);
         count++;
         GenBank.GBIterator++;
     }
    for(int i=0;i<SampleNum;i++)
    {
        dist_calc(i);
        EM(i);
        for(int j=0;j<TaxonomyID.size();j++)
        {
            if(InitialGuess[j]!=0)
            {
                map<string, double*>::iterator it=Abundance.find(rTaxonomyID[j]);
                if(it==Abundance.end())
                {
                    Abundance[rTaxonomyID[j]]=new double[SampleNum];
                    memset(Abundance[rTaxonomyID[j]],0,sizeof(double)*SampleNum);
                    Abundance[rTaxonomyID[j]][i]=InitialGuess[j];
                }
                else
                {
                    Abundance[rTaxonomyID[j]][i]=InitialGuess[j];
                }
            }
        }
        delete[] InitialGuess;
    }
    return 1;
}

int AnalysisPipeline::open(string* BLASTFileNames, string GenBankFileName, bool BuildTaxonFlag, int OperationFlag, bool EMFlag, bool WilcoxonFlag, bool MedianNormalizationFlag)
{
    if(BuildTaxonFlag==true)
    {
        if(!GenBank.extract_taxonomy(GenBankFileName+".gen"))
            return 0;
        GenBank.output_taxonomy(GenBankFileName+".gen");
        GenBank.load_taxonomy(GenBankFileName+".gen");
    }
    else if(!GenBank.load_taxonomy(GenBankFileName+".gen"))
    {
        if(!GenBank.extract_taxonomy(GenBankFileName+".gen"))
            return 0;
        GenBank.output_taxonomy(GenBankFileName+".gen");
        GenBank.load_taxonomy(GenBankFileName+".gen");
    }
    RDPQuery.load(GenBankFileName+".fa");   
    BLASTMapping=new BLASTResult[SampleNum];
    for(int i=0;i<SampleNum;i++)
    {
        GenBank.Iterator_reset(OperationFlag);
        while(!GenBank.Iterator_end(OperationFlag))
        {
            BLASTMapping[i].DBAbundance_set((*GenBank.GBIterator).first,(double)(*GenBank.GBIterator).second/((1-(*GenBank.GBIterator).second)*MIN(OperationFlag,1)+(*GenBank.GBIterator).second));  
            GenBank.GBIterator++;
        }
        BLASTMapping[i].SetThresholds(BLASTEvalueThreshold, BLASTIdentityThreshold, BLASTMappingLengthThreshold);
        BLASTMapping[i].Import(BLASTFileNames[i],GenBank, RDPQuery, OperationFlag, TaxonomicLevel);
    }
    if(EMFlag==true)
        infer_taxonomy_EM();
    else
        infer_taxonomy(OperationFlag, MedianNormalizationFlag);
    statistics(WilcoxonFlag);
    return 1;    
}

int AnalysisPipeline::output(string OutputFileName, bool WilcoxonFlag)
{
    fstream OutputFile, StatOutputFile;
    OutputFile.open(OutputFileName.c_str(),ios_base::out);
    if(OutputFile.fail()){cout<<"cannot open output file"; return 0;}
    for(int i=0;i<SampleNum;i++)
    {
        OutputFile<<AmbiSize[i];
        if(i==SampleNum-1)
            OutputFile<<"\n";
        else
            OutputFile<<"\t";
    }
    map<string,double*>::iterator it=Abundance.begin();
    while(it!=Abundance.end())
    {
        OutputFile<<(*it).first<<"\t";
        for(int i=0;i<SampleNum;i++)
        {
            OutputFile<<(*it).second[i];
            if(i==SampleNum-1)
                OutputFile<<"\n";
            else
                OutputFile<<"\t";
        }
         it++;
    }
    OutputFile.close();

//*********************************************************************************
    StatOutputFile.open((OutputFileName+".wilcoxon").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open wilcoxon test output file"; return 0;}
    it=WilcoxonScore.begin();
    while(it!=WilcoxonScore.end())
    {
        StatOutputFile<<(*it).first<<"\t";
        for(int j=0;j<6;j++)
        {
            StatOutputFile<<(*it).second[j];
            if(j==5)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
        it++;
    }
    StatOutputFile.close();
    
//********************************************************************************************  

//*********************************************************************************
    if(WilcoxonFlag)
    {
    StatOutputFile.open((OutputFileName+".wilcoxon_special_sensitivity").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open wilcoxon rank sum test output file"; return 0;}
    for(int i=0;i<6;i++)
    {
        for(int j=0;j<9;j++)
        {
            for(int k=0;k<5;k++)
            {
                 if(WilcoxonMatrixSensitivity[k][i][j]==1.0)
                 {
                     StatOutputFile<<k;
                     if(k==0)
                         StatOutputFile<<"-";
                     if(k!=0&&WilcoxonMatrixSensitivity[k-1][i][j]==0.0)
                         StatOutputFile<<"=";
                     if(k!=0&&WilcoxonMatrixSensitivity[k-1][i][j]==-1.0)
                         StatOutputFile<<"+";
                     break;
                 }
                 if(k==4)
                 {
                     StatOutputFile<<k+1;
                     if(WilcoxonMatrixSensitivity[k][i][j]==0.0)
                         StatOutputFile<<"=";
                     else
                         StatOutputFile<<"+";
                 }
            }
            if(j==8)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
    }         
    StatOutputFile.close();
    for(int i=0;i<5;i++)
    {
       for(int j=0;j<6;j++)
           delete[] WilcoxonMatrixSensitivity[i][j];
       delete[] WilcoxonMatrixSensitivity[i];
    }
    delete[] WilcoxonMatrixSensitivity;

    StatOutputFile.open((OutputFileName+".wilcoxon_special_specificity").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open wilcoxon rank sum test output file"; return 0;}
    for(int i=5;i>=0;i--)
    {
        for(int j=8;j>=0;j--)
        {
            for(int k=4;k>=0;k--)
            {
                 if(WilcoxonMatrixSpecificity[k][i][j]==1.0)
                 {
                     StatOutputFile<<k+1;
                     if(k==4)
                         StatOutputFile<<"-";
                     if(k!=4&&WilcoxonMatrixSpecificity[k+1][i][j]==0.0)
                         StatOutputFile<<"=";
                     if(k!=4&&WilcoxonMatrixSpecificity[k+1][i][j]==-1.0)
                         StatOutputFile<<"+";                     
                     break;
                 }
                 if(k==0)
                 {
                     StatOutputFile<<k;
                     if(WilcoxonMatrixSpecificity[k][i][j]==0.0)
                         StatOutputFile<<"=";
                     else
                         StatOutputFile<<"+";
                 }
            }
            if(j==0)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
    }         
    StatOutputFile.close();
    for(int i=0;i<5;i++)
    {
       for(int j=0;j<6;j++)
           delete[] WilcoxonMatrixSpecificity[i][j];
       delete[] WilcoxonMatrixSpecificity[i];
    }
    delete[] WilcoxonMatrixSpecificity;

    StatOutputFile.open((OutputFileName+".wilcoxon_special_bray_curtis").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open wilcoxon rank sum test output file"; return 0;}
    for(int i=0;i<6;i++)
    {
        for(int j=0;j<9;j++)
        {
            for(int k=0;k<5;k++)
            {
                 if(WilcoxonMatrixBrayCurtis[k][i][j]==-1.0)
                 {
                     StatOutputFile<<k;
                     if(k==0)
                         StatOutputFile<<"-";
                     if(k!=0&&WilcoxonMatrixBrayCurtis[k-1][i][j]==0.0)
                         StatOutputFile<<"=";
                     if(k!=0&&WilcoxonMatrixBrayCurtis[k-1][i][j]==1.0)
                         StatOutputFile<<"+";                         
                     break;
                 }
                 if(k==4)
                 {
                     StatOutputFile<<k+1;
                     if(WilcoxonMatrixBrayCurtis[k][i][j]==0.0)
                         StatOutputFile<<"=";
                     else
                         StatOutputFile<<"+";
                 }
            }
            if(j==8)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
    }         
    StatOutputFile.close();
    for(int i=0;i<5;i++)
    {
       for(int j=0;j<6;j++)
           delete[] WilcoxonMatrixBrayCurtis[i][j];
       delete[] WilcoxonMatrixBrayCurtis[i];
    }
    delete[] WilcoxonMatrixBrayCurtis;
    }    // if WilcoxonFlag
//********************************************************************************************  
    
//***************************************************************************************    
    StatOutputFile.open((OutputFileName+".jaccard").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open jaccard index output file"; return 0;}
    for(int i=0;i<SampleNum;i++)
        for(int j=0;j<SampleNum;j++)
        {
            StatOutputFile<<JaccardIndexValue[i][j];
            if(j==SampleNum-1)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
    StatOutputFile.close();

//*****************************************************************************************
    StatOutputFile.open((OutputFileName+".jaccard_special").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open jaccard index output file"; return 0;}
    for(int i=1;i<SampleNum;i++)
    {
        StatOutputFile<<JaccardIndexValue[0][i];
        if(i%10==0)
            StatOutputFile<<"\n";
        else
            StatOutputFile<<"\t";
        //if((i-1)%120==0)
            //StatOutputFile<<"\n";
    }
    StatOutputFile.close();
//*****************************************************************************************
    
    StatOutputFile.open((OutputFileName+".sensitivity").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open sensitivity output file"; return 0;}
    for(int i=0;i<SampleNum;i++)
        for(int j=0;j<SampleNum;j++)
        {
            StatOutputFile<<SensitivityMatrix[i][j];
            if(j==SampleNum-1)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
    StatOutputFile.close();    

//*****************************************************************************************
    StatOutputFile.open((OutputFileName+".sensitivity_special").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open sensitivity output file"; return 0;}
    for(int i=1;i<SampleNum;i++)
    {
        StatOutputFile<<SensitivityMatrix[0][i];
        if(i%10==0)
            StatOutputFile<<"\n";
        else
            StatOutputFile<<"\t";
        //if((i-1)%120==0)
            //StatOutputFile<<"\n";
    }
    for(int i=1;i<SampleNum;i++)
    {
        StatOutputFile<<SensitivityMatrix[i][0];
        if(i%10==0)
            StatOutputFile<<"\n";
        else
            StatOutputFile<<"\t";
        //if((i-1)%120==0)
            //StatOutputFile<<"\n";
    }
    StatOutputFile.close();
//*******************************************************************************************

    StatOutputFile.open((OutputFileName+".theta").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open theta index output file"; return 0;}
    for(int i=0;i<SampleNum;i++)
        for(int j=0;j<SampleNum;j++)
        {
            StatOutputFile<<ThetaIndexValue[i][j];
            if(j==SampleNum-1)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
    StatOutputFile.close();

//********************************************************************************************
    StatOutputFile.open((OutputFileName+".theta_special").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open theta index output file"; return 0;}
    for(int i=1;i<SampleNum;i++)
    {
        StatOutputFile<<ThetaIndexValue[0][i];
        if(i%10==0)
            StatOutputFile<<"\n";
        else
            StatOutputFile<<"\t";
        //if((i-1)%120==0)
            //StatOutputFile<<"\n";
    }
    StatOutputFile.close();
//********************************************************************************************   
 
    StatOutputFile.open((OutputFileName+".bray_curtis").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open theta index output file"; return 0;}
    for(int i=0;i<SampleNum;i++)
        for(int j=0;j<SampleNum;j++)
        {
            StatOutputFile<<BrayCurtisDissimilarity[i][j];
            if(j==SampleNum-1)
                StatOutputFile<<"\n";
            else
                StatOutputFile<<"\t";
        }
    StatOutputFile.close();
 
//********************************************************************************************
    StatOutputFile.open((OutputFileName+".bray_curtis_special").c_str(),ios_base::out);
    if(StatOutputFile.fail()){cout<<"cannot open theta index output file"; return 0;}
    for(int i=1;i<SampleNum;i++)
    {
        StatOutputFile<<BrayCurtisDissimilarity[0][i];
        if(i%10==0)
            StatOutputFile<<"\n";
        else
            StatOutputFile<<"\t";
        //if((i-1)%120==0)
            //StatOutputFile<<"\n";
    }
    StatOutputFile.close();
//********************************************************************************************  

    return 1;
}

int  AnalysisPipeline::statistics(bool WilcoxonFlag)
{
    SensitivityMatrix=new double*[SampleNum];
    JaccardIndexValue=new double*[SampleNum];
    ThetaIndexValue=new double*[SampleNum];
    BrayCurtisDissimilarity=new double*[SampleNum];
    for(int i=0;i<SampleNum;i++)
    {
        SensitivityMatrix[i]=new double[SampleNum];
        JaccardIndexValue[i]=new double[SampleNum];
        memset(JaccardIndexValue[i],0,sizeof(double)*SampleNum);        
        ThetaIndexValue[i]=new double[SampleNum];
        memset(ThetaIndexValue[i],0,sizeof(double)*SampleNum);
        BrayCurtisDissimilarity[i]=new double[SampleNum];
        memset(BrayCurtisDissimilarity[i],0,sizeof(double)*SampleNum);        
    }
    SensitivityMatrix=Sensitive.calc(Abundance, SampleNum, SensitivityMatrix);
    cout<<"Sensitivity done\n";
    JaccardIndexValue=Jaccard.calc(Abundance, SampleNum, JaccardIndexValue);
    cout<<"Jaccard Index done\n";
    ThetaIndexValue=Theta.calc(Abundance, SampleNum, ThetaIndexValue);
    cout<<"Theta Index done\n";
    BrayCurtisDissimilarity=BCurtis.calc(Abundance, SampleNum, BrayCurtisDissimilarity);
    cout<<"Bray-Curtis done\n";
    
    if(WilcoxonFlag)
    {
    //***************************** Wilcoxon for Simulated Data ***************************
    //Sensitivity
    double temparray1[10], temparray2[10];
    WilcoxonMatrixSensitivity=new double**[5];
    for(int i=0;i<5;i++)
    {
        WilcoxonMatrixSensitivity[i]=new double*[6];
        for(int t=0;t<10;t++)
            temparray1[t]=SensitivityMatrix[0][10*i+t+1];
        for(int j=0;j<6;j++)
        {
            WilcoxonMatrixSensitivity[i][j]=new double[9];
            for(int k=0;k<9;k++)
            {
                for(int t=0;t<10;t++)
                    temparray2[t]=SensitivityMatrix[0][10*(9*j+k+5)+1+t];
                WilcoxonMatrixSensitivity[i][j][k]=Wil.wilcoxon_rank_sum_test(temparray1, temparray2, 10, 10);
                if(WilcoxonMatrixSensitivity[i][j][k]!=0.0&&WilcoxonMatrixSensitivity[i][j][k]!=1.0&&WilcoxonMatrixSensitivity[i][j][k]!=-1.0)
                    cout<<"Sensitivity\t"<<i<<"\t"<<j<<"\t"<<k<<"\tis not 0, -1 or 1\n";
            }
        }
    }
    //Specificity
    WilcoxonMatrixSpecificity=new double**[5];
    for(int i=0;i<5;i++)
    {
        WilcoxonMatrixSpecificity[i]=new double*[6];
        for(int t=0;t<10;t++)
            temparray1[t]=SensitivityMatrix[10*i+t+1][0];
        for(int j=0;j<6;j++)
        {
            WilcoxonMatrixSpecificity[i][j]=new double[9];
            for(int k=0;k<9;k++)
            {
                for(int t=0;t<10;t++)
                    temparray2[t]=SensitivityMatrix[10*(9*j+k+5)+1+t][0];
                WilcoxonMatrixSpecificity[i][j][k]=Wil.wilcoxon_rank_sum_test(temparray1, temparray2, 10, 10);
                if(WilcoxonMatrixSpecificity[i][j][k]!=0.0&&WilcoxonMatrixSpecificity[i][j][k]!=1.0&&WilcoxonMatrixSpecificity[i][j][k]!=-1.0)
                    cout<<"Specificity\t"<<i<<"\t"<<j<<"\t"<<k<<"\tis not 0, -1 or 1\n";
            }
        }
    }
    //Bray-Curtis
    WilcoxonMatrixBrayCurtis=new double**[5];
    for(int i=0;i<5;i++)
    {
        WilcoxonMatrixBrayCurtis[i]=new double*[6];
        for(int t=0;t<10;t++)
            temparray1[t]=BrayCurtisDissimilarity[0][10*i+t+1];
        for(int j=0;j<6;j++)
        {
            WilcoxonMatrixBrayCurtis[i][j]=new double[9];
            for(int k=0;k<9;k++)
            {
                for(int t=0;t<10;t++)
                    temparray2[t]=BrayCurtisDissimilarity[0][10*(9*j+k+5)+1+t];
                WilcoxonMatrixBrayCurtis[i][j][k]=Wil.wilcoxon_rank_sum_test(temparray1, temparray2, 10, 10);
                if(WilcoxonMatrixBrayCurtis[i][j][k]!=0.0&&WilcoxonMatrixBrayCurtis[i][j][k]!=1.0&&WilcoxonMatrixBrayCurtis[i][j][k]!=-1.0)
                    cout<<"BrayCurtis\t"<<i<<"\t"<<j<<"\t"<<k<<"\tis not 0, -1 or 1\n";
            }
        }
    }            
    /*
    //***************************** Wilcoxon Tests *********************************
    double** temparray=new double*[4];
    for(int i=0;i<4;i++)
        temparray[i]=new double[18];
    map<string, double*>::iterator it=Abundance.begin();
    while(it!=Abundance.end())
    {
        WilcoxonScore[(*it).first]=new double[6];
        for(int i=0;i<4;i++)
            for(int j=0;j<18;j++)
                temparray[i][j]=(*it).second[j+18*i];
        int count=0;
        for(int i=0;i<4;i++)
            for(int j=i+1;j<4;j++)
            {
               WilcoxonScore[(*it).first][count]=Wil.wilcoxon_signed_rank_test(temparray[i], temparray[j], 18);
               count++;
            }
        it++;
    }
    for(int i=0;i<4;i++)
        delete[] temparray[i];
    delete[] temparray;
    //**********************************************************************************
    */
    }
    cout<<"Wilcoxon done\n";
    return 1;
}

void  AnalysisPipeline::normalize_median()
{
    double** tVector=new double*[SampleNum];
    for(int i=0;i<SampleNum;i++)
    {
        tVector[i]=new double[Abundance.size()];
        memset(tVector[i], 0, sizeof(double)*Abundance.size());
    }
    Median=new double[SampleNum];
    memset(Median, 0, sizeof(double)*SampleNum);
    map<string, double*>::iterator it=Abundance.begin();
    int count=0;
    while(it!=Abundance.end())
    {
        for(int i=0;i<SampleNum;i++)
            tVector[i][count]=(*it).second[i];
        it++;
        count++;
    }
    for(int i=0;i<SampleNum;i++)
        Median[i]=find_median(tVector[i], Abundance.size());
    double* tempMedian=new double[SampleNum];
    memcpy(tempMedian, Median, sizeof(double)*SampleNum);
    double tMedian=find_median(tempMedian, SampleNum);
    it=Abundance.begin();
    while(it!=Abundance.end())
    {
        for(int i=0;i<SampleNum;i++)
            (*it).second[i]=(double)(*it).second[i]*tMedian/Median[i];
        it++;
    }
    for(int i=0;i<SampleNum;i++)
        cout<<Median[i]<<"\t";
    cout<<"\n";
    for(int i=0;i<SampleNum;i++)
        delete[] tVector[i];
    delete[] tVector;
    delete[] tempMedian;
    return;
}

double AnalysisPipeline::find_median(double* vector, int size)
{
    for(int i=0;i<size;i++)
    {
        for(int j=i+1;j<size;j++)
            if(vector[i]>vector[j])
                Wil.swap(vector[i],vector[j]);
    }
    int count=0;
    for(int i=0;i<size;i++)
        if(vector[i]==0)
            count++;
    if((size-count)%2==0)
            return (vector[count+(size-count)/2]+vector[count+(size-count)/2-1])/2;
    if((size-count)%2!=0)
            return vector[count+(size-count-1)/2];
    return 0;
}

int AnalysisPipeline::infer_taxonomy_new()
{    
}
